using DataFrames, Plots, Statistics, JSON, Clp
plotlyjs()
include(joinpath(dirname(pwd()),"src/TuLiPa.jl")); # using Dates, JuMP, HiGHS, CSV, Clustering
include("Demo 1a - Deterministic power market with dummy data.jl"); # Hid some functions here to make the notebook more readable
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
┌ Warning: backend `PlotlyJS` is not installed. └ @ Plots C:\Users\jgrc\.julia\packages\Plots\bMtsB\src\backends.jl:43
We show how you can set up a simple deterministic power market problem with TuLiPa. The problem will minimize the cost of covering the demand with production from different technologies.
The problem contains:
These elements need their own model objects. Each model object (e.g. hydropower production or the water balance) has an id, horizon (with time resolution) and special properties. The properties of the object contain information about which elements must be built and updated in the LP problem, as well as a dataset. The dataset can contain values for different weather scenarios, and different model years (e.g. installed capacity in 2021 or 2025). When we update the problem with a weather scenario (scenario time) and a model year (data time), the model objects will build the LP problem for this input data. Then the problem can be solved.
In this demo we use dummy data, and profiles for wind, solar, demand and inflow from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/
The model objects are assembled from a list of dataelements. The dataelements have a concept, concepttype, unique instancename and a data field. The conceptname and concepttype decide how the dataelement should be converted to a model object, while the data field holds all the data that should be put into the model object.
The data elements are stored in a flat structure where each data element holds its information and possibly references to other dataelements they are linked to (for example can an upper capacity hold references to its parameter (stored in another data element), and references to the variable this upper capacity should apply to).
# We have hid these functions in the file Demo 1a - Deterministic power market with dummy data.jl
# Alternatively take a look at Demo 1b where all the code is included in the notebook
elements = gettestdataset();
println(string("Number of elements: ", length(elements)))
printdicts(elements, 10) # print first 10 elements
# printdicts(elements) # print all elements
Number of elements: 108
[
{
"conceptname": "Balance",
"typename": "BaseBalance",
"instancename": "PowerBalance_NO2",
"value": {
"Commodity": "Power"
}
},
{
"conceptname": "Flow",
"typename": "BaseFlow",
"instancename": "SlackVarPowerBalance_NO2",
"value": {}
},
{
"conceptname": "Arrow",
"typename": "BaseArrow",
"instancename": "SlackArrowPowerBalance_NO2",
"value": {
"Balance": "PowerBalance_NO2",
"Flow": "SlackVarPowerBalance_NO2",
"Conversion": 1.0,
"Direction": "Out"
}
},
{
"conceptname": "Balance",
"typename": "BaseBalance",
"instancename": "PowerBalance_GER",
"value": {
"Commodity": "Power"
}
},
{
"conceptname": "Flow",
"typename": "BaseFlow",
"instancename": "SlackVarPowerBalance_GER",
"value": {}
},
{
"conceptname": "Arrow",
"typename": "BaseArrow",
"instancename": "SlackArrowPowerBalance_GER",
"value": {
"Balance": "PowerBalance_GER",
"Flow": "SlackVarPowerBalance_GER",
"Conversion": 1.0,
"Direction": "Out"
}
},
{
"conceptname": "Flow",
"typename": "BaseFlow",
"instancename": "PowerBalance_NO2->PowerBalance_GER",
"value": {}
},
{
"conceptname": "Arrow",
"typename": "BaseArrow",
"instancename": "PowerBalance_NO2->PowerBalance_GERFrom",
"value": {
"Balance": "PowerBalance_NO2",
"Flow": "PowerBalance_NO2->PowerBalance_GER",
"Conversion": 1.0,
"Direction": "Out"
}
},
{
"conceptname": "Arrow",
"typename": "BaseArrow",
"instancename": "PowerBalance_NO2->PowerBalance_GERTo",
"value": {
"Balance": "PowerBalance_GER",
"Flow": "PowerBalance_NO2->PowerBalance_GER",
"Conversion": 0.97,
"Direction": "In"
}
},
{
"conceptname": "Capacity",
"typename": "PositiveCapacity",
"instancename": "CapPowerBalance_NO2->PowerBalance_GER",
"value": {
"Param": "CapPowerBalance_NO2->PowerBalance_GER",
"WhichInstance": "PowerBalance_NO2->PowerBalance_GER",
"WhichConcept": "Flow",
"Bound": "Upper"
}
}
]
Define some model settings that completes the dataset.
# Select the weather scenarios 1981-1982 (some of the profiles are 1981-2010, we only want to use two years)
# Used by RotatingTimeVector to reduce the amount of data
scenarioyearstart = 1981
scenarioyearstop = 1983
addscenariotimeperiod!(elements, "ScenarioTimePeriod", getisoyearstart(scenarioyearstart), getisoyearstart(scenarioyearstop));
# We choose the horizon (time-resolution) of the commodities. We set the duration of the horizons to 3 years
# We want the variables connected to power (daily) to be more detailed than the hydro variables (weekly)
power_horizon = SequentialHorizon(364*3, Day(1))
hydro_horizon = SequentialHorizon(52*3, Week(1))
set_horizon!(elements, "Power", power_horizon)
set_horizon!(elements, "Hydro", hydro_horizon);
# Storages have state-dependant variables that need a boundary condition
# We set the starting storage to be equal to the ending storage, x[0] = x[T] (for horizon where t in 1:T)
push!(elements, getelement(BOUNDARYCONDITION_CONCEPT, "StartEqualStop", "StartEqualStop_StorageResNO2",
(WHICHINSTANCE, "StorageResNO2"),
(WHICHCONCEPT, STORAGE_CONCEPT)));
modelobjects = getmodelobjects(elements)
keys(modelobjects)
KeySet for a Dict{Any, Any} with 15 entries. Keys:
Id("Flow", "SlackVarPowerBalance_NO2")
Id("Flow", "GenCoalGER")
Id("Balance", "PowerBalance_GER")
Id("Flow", "ReleaseResNO2")
Id("Flow", "PowerBalance_GER->PowerBalance_NO2")
Id("Balance", "ResNO2")
Id("Flow", "GenBioGER")
Id("Flow", "GenNucGER")
Id("Balance", "PowerBalance_NO2")
Id("Flow", "SpillResNO2")
Id("Flow", "SlackVarPowerBalance_GER")
Id("Flow", "GenGasGER")
Id("BoundaryCondition", "StartEqualStop_StorageResNO2")
Id("Storage", "StorageResNO2")
Id("Flow", "PowerBalance_NO2->PowerBalance_GER")
println("We display the fields in the struct")
display(fieldnames(typeof(modelobjects[Id("Flow", "ReleaseResNO2")])))
println("\nThe Id shows that the struct is of type Flow")
display(getid(modelobjects[Id("Flow", "ReleaseResNO2")]))
println("\nThe Horizon consist of 1092 periods with timedelta of 1 day (total 3 years)")
display(gethorizon(modelobjects[Id("Flow", "ReleaseResNO2")]))
println("\nThe lower bound is LowerZeroCapacity so that the variable is non-negative")
display(getlb(modelobjects[Id("Flow", "ReleaseResNO2")]))
println("\nThe upper bound is 6000 m3/s, which is stored in a M3SToMM3SeriesParam. M3SToMM3SeriesParam converts the upper bound into Mm3 depending on the horizon of the problem")
display(getub(modelobjects[Id("Flow", "ReleaseResNO2")]))
println("\nThere are two arrows. The first says that the release takes water out from the water balance and the second says that it contributes into the power balance")
for arrow in getarrows(modelobjects[Id("Flow", "ReleaseResNO2")])
display(getid(arrow))
# display(arrow)
end
We display the fields in the struct
(:id, :horizon, :ub, :lb, :costs, :sumcost, :arrows, :metadata)
The Id shows that the struct is of type Flow
Id("Flow", "ReleaseResNO2")
The Horizon consist of 1092 periods with timedelta of 1 day (total 3 years)
SequentialHorizon(SequentialPeriods([(1092, Millisecond(86400000))]), nothing)
The lower bound is LowerZeroCapacity so that the variable is non-negative
LowerZeroCapacity()
The upper bound is 6000 m3/s, which is stored in a M3SToMM3SeriesParam. M3SToMM3SeriesParam converts the upper bound into Mm3 depending on the horizon of the problem
PositiveCapacity(Id("Capacity", "CapReleaseResNO2"), M3SToMM3SeriesParam{ConstantTimeVector, ConstantTimeVector}(ConstantTimeVector(6000.0), ConstantTimeVector(1.0)), true)
There are two arrows. The first says that the release takes water out from the water balance and the second says that it contributes into the power balance
Id("Arrow", "WaterDemandReleaseResNO2")
Id("Arrow", "PowerSupplyReleaseResNO2")
# It is stored in one of the rhsterms of PowerBalance_NO2
rhsterm = getrhsterms(modelobjects[Id("Balance", "PowerBalance_GER")])[1]
# We start with the installed capacity in 2023 and the weather scenario 1981
start = TwoTime(getisoyearstart(2023), getisoyearstart(1981))
# We look at the weekly production in GWh for three years from the starting point
horizon = SequentialHorizon(52*3, Week(1))
# Vectors to store results
i = Vector{typeof(getdatatime(start))}(undef, getnumperiods(horizon))
v = zeros(getnumperiods(horizon))
# Query results for each period in the horizon
for t in 1:getnumperiods(horizon)
querystart = getstarttime(horizon, t, start)
querydelta = gettimedelta(horizon, t)
i[t] = getdatatime(querystart)
v[t] = getparamvalue(rhsterm, querystart, querydelta)
end
# Plot
display(plot(i, v, legend=false, size=(700,300), ylabel = "GWh", title = "Wind power production GER"))
We build a problem object for a minimization problem. The problem object consist of:
The problem is also initialized. Variables and balances are built, and variables are put into balances and the objective function if coefficients are constant regardless of scenario or horizon period.
mymodel = JuMP.Model(Clp.Optimizer)
set_silent(mymodel)
prob = JuMP_Prob(modelobjects, mymodel)
# prob = HiGHS_Prob(modelobjects) # we could also have used HiGHS_Prob, but printing the mathematical model is not supported there
prob.model
A JuMP Model
Minimization problem with:
Variables: 10141
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 2341 constraints
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 10140 constraints
`VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 7800 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Clp
Names registered in the model: BalancePowerBalance_GER, BalancePowerBalance_NO2, BalanceResNO2, BoundaryConditionEqStartEqualStop_StorageResNO2, FlowGenBioGER, FlowGenCoalGER, FlowGenGasGER, FlowGenNucGER, FlowPowerBalance_GER->PowerBalance_NO2, FlowPowerBalance_NO2->PowerBalance_GER, FlowReleaseResNO2, FlowSlackVarPowerBalance_GER, FlowSlackVarPowerBalance_NO2, FlowSpillResNO2, StorageStartStorageResNO2, StorageStorageResNO2
The hydrorbalance has been initialized. But the inflow is not added since the value depends on the chosen scenario and period in the horizon.
display(prob.model[:BalanceResNO2][1])
Variables are now put into balances and the objective function if coefficients are dependant on the scenario or horizon period. RHSTerms like the inflow are also updated, but it is stored in the problem object rather than updated in the JuMP model straight away. When solve! is called, the problem object will check if any of the RHSTerm has been updated. If that is the case, the JuMP model will be updated with the RHSTerms.
t = TwoTime(getisoyearstart(2021), getisoyearstart(1981))
update!(prob, t)
@time solve!(prob)
getobjectivevalue(prob)
5.450143 seconds (9.89 M allocations: 553.758 MiB, 5.17% gc time, 98.80% compilation time)
9.61646303403377e10
Now the hydropower balance also includes the inflow, which is 156 Mm3 in the second week of the problem.
display(prob.model[:BalanceResNO2][1])
t = TwoTime(getisoyearstart(2024), getisoyearstart(1982))
update!(prob, t)
solve!(prob)
display(getobjectivevalue(prob))
display(prob.model[:BalanceResNO2][1])
8.194441022906743e10
The prices are in EUR/MWh. The price decline from 2025 is a result of the new gas price level in 2025.
plot(-dual.(prob.model[:BalancePowerBalance_GER])/1000)
# plot(value.(prob.model[:FlowGenGasGER])) # Or we could get the daily gas production in GWh
numperiods = getnumperiods(power_horizon)
x = [getdatatime(t) + getstartduration(power_horizon, s) for s in 1:numperiods]
y = [-getcondual(prob, Id("Balance","PowerBalance_GER"), s) for s in 1:numperiods]/1000
plot(x,y)
The function runs several scenarios and plot prices, supplies, demands and hydro storages for each scenario. The function only plots results for chosen price areas.
In Demo 2 a more complex version of this function is used, which support more types of objects. This includes AdaptiveHorizons (horizon built based on residual load), Horizons with offset, SegmentedArrow (for hydropower plants with PQ-curves), BaseAggSupplyCurve (for aggregated power plants) and ExogenBalance (exogen price areas)
# We have hid this function in the julia file: Demo 1a - Deterministic power market with dummy data.jl
We choose scenarios that start in the model year 2024 and weather scenarios that start in 1981/1982/1983
We get expected results for the deterministic market problem
scenarios = getscenarios(2024; years=1981:1983)
runscenarios(scenarios, modelobjects)
# alternative for only seeing results from one price area
# resultobjects = getpowerobjects(modelobjects,["NO2"])
# runscenarios(scenarios, modelobjects, resultobjects)
Objective value in scenario 1: 8.264928001510184e10 Objective value in scenario 2: 8.19444102290675e10 Objective value in scenario 3: 8.303707657019687e10
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 131.058 DemandGER 468.054 2 │ WindGER 111.944 DemandNO2 48.4256 3 │ SolarGER 70.6705 PowerBalance_NO2->PowerBalance_G… 8.62053 4 │ GenCoalGER 60.4118 SlackVarPowerBalance_NO2 1.80325 5 │ GenNucGER 43.5204 PowerBalance_GER->PowerBalance_N… 0.984684 6 │ GenBioGER 43.1119 SlackVarPowerBalance_GER 0.0400589 7 │ ReleaseResNO2 36.954 missing missing 8 │ RoRNO2 15.9198 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.36191 missing missing 10 │ WindNO2 5.02039 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.955143 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 527.929 527.929
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 129.634 DemandGER 467.293 2 │ WindGER 109.393 DemandNO2 47.626 3 │ SolarGER 73.102 PowerBalance_NO2->PowerBalance_G… 9.10136 4 │ GenCoalGER 60.1686 SlackVarPowerBalance_NO2 1.48872 5 │ GenNucGER 43.5209 PowerBalance_GER->PowerBalance_N… 0.412345 6 │ GenBioGER 43.1053 SlackVarPowerBalance_GER 0.0463198 7 │ ReleaseResNO2 37.2763 missing missing 8 │ RoRNO2 15.4017 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.82832 missing missing 10 │ WindNO2 5.13809 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.399974 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 525.968 525.968
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 131.754 DemandGER 468.242 2 │ WindGER 111.142 DemandNO2 48.4851 3 │ SolarGER 70.6877 PowerBalance_NO2->PowerBalance_G… 8.56625 4 │ GenCoalGER 60.7569 SlackVarPowerBalance_NO2 1.79717 5 │ GenNucGER 43.5204 PowerBalance_GER->PowerBalance_N… 0.994443 6 │ GenBioGER 43.1119 SlackVarPowerBalance_GER 0.0463198 7 │ ReleaseResNO2 36.9577 missing missing 8 │ RoRNO2 15.9212 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.30926 missing missing 10 │ WindNO2 5.005 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.96461 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 528.131 528.131